---
title: "R2.TxTrends"
output: pdf_document
---

# Load Packages
```{r}
library(tidyverse)
library(estimatr)
```

# Load Data
```{r}
load("Data/data_final.Rdata")
```

# Set Covariates
```{r}
covs_list <- c("tx_Hi", "as.factor(DistN)", "ln_PopDens_70", "Cities_10km", "ln_Confs_All", "PrRecH_Adm3", "Isolados1950", "PrWheatAr_1972", "LandIneqRat", "TerrainSlopeIndex", "MedAltitude", "Communist_1975", "MinDistPCP_Map")
```

# Regress Land Expropriation ~ Treatment + Covariates after changing the treatment threshold from the 10th to 90th percentile (in increments of 5) of the mismatch index (resid). 
```{r}

thresholds <- seq(0.1, 0.9, by = 0.05) # percentiles
holder <- list()
for (i in 1:length(thresholds)){ # i<-1
data_subset_tx <- data_final %>% select(-tx_Hi) %>% ungroup() %>%
  mutate(tx_Hi = ifelse(resid >= quantile(resid, probs = thresholds[i], na.rm=T), 1, 0)) %>%
  filter(!is.na(tx_Hi))

holder[[i]]<- lm_robust(as.formula(paste("pr_ha_ord", paste(covs_list, collapse=" + "), sep = " ~ ")), 
               data = data_subset_tx)

}

output <- data.frame(thresholds, tx_Hi = NA, se = NA, n = NA)

for (l in 1:length(holder)){
  output[l,"tx_Hi"]<-holder[[l]]$coefficients["tx_Hi"]
  output[l,"se"]<-holder[[l]]$std.error["tx_Hi"]
  output[l,"n"]<-length(holder[[l]]$fitted.values)
}

ggplot(output, aes(x = thresholds)) +
  geom_hline(aes(yintercept = 0), col = "red") +  xlim(0.1, 0.9) +
  geom_point(aes(y = tx_Hi)) + 
  geom_linerange(aes(ymin=tx_Hi - (qt(0.975, n)*se), ymax=tx_Hi + (qt(0.975, n)*se)), size = 0.5) +
  geom_linerange(aes(ymin=tx_Hi - (qt(0.95, n)*se), ymax=tx_Hi + (qt(0.95, n)*se)), size = 0.75) +
  scale_x_continuous(breaks = seq(0.1:0.9, by = 0.1)) +
  xlab("Degree of Mismatch (Percentile), Dichotomized") + ylab("Effect of Mismatch on % Hectares Expropriated") +
  theme_light()

ggsave(plot = last_plot(),
       "Figures/fg3a_TxExprop.pdf", 
       width = 5, height = 5)

```